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Abstract. - A transfer matrix approach to study ballistic charge transport in few-layer graphene 
with chiral-symmetric stacking configurations is developed. We demonstrate that the chiral sym- 
metry justifies a non-Abelian gauge transformation at the spectral degeneracy point (zero energy). 
This transformation proves the equivalence of zero-energy transport properties of the multilayer 
to those of the system of uncoupled monolayers. Similar transformation can be applied in order 
to gauge away an arbitrary magnetic field, weak strain, and hopping disorder in the bulk of the 
sample. Finally, we calculate the full-counting statistics at arbitrary energy for different stacking 
configurations. The predicted gate-voltage dependence of conductance and noise can be measured 
in clean multilayer samples with generic metallic leads. 



^ ' Introduction. — Discovery of graphene [T] has revived 
Q an interest to a wide range of phenomena that are ex- 
Q plicitly or implicitly linked to the symmetry of the band 
'~~btructure [2l|3]. Graphene and topological metals realised 
in BiSb compounds [3] are very recent examples where ex- 
^ otic band-structure properties are responsible for unusual 
quantum effects in the charge transport. The well-known 
one is the integer quantum Hall effect that takes on re- 
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iharkably distinct forms in mono- and bi-layer graphene 



.[SHZ]- The continuing development of device fabrication 
technology promises interesting applications of these ef- 
Q^fects in microelectronics, chemical sensing, and quantum 
information. 

IL" ' At the same time there is a growing interest in new 
.^functionalities of few-layer graphene that is inspired by 
^the proposals to control valley polarisation [HI IS], band 
Hgaps [TU], and magnetic exchange [TT] with the help of 



gate electrodes. The experimental evidence supporting 
these theoretical proposals is currently building up [T2IIT3] . 
Few-layer graphene is also a prime candidate for the obser- 
vation of excitonic condensation caused by an attractive 
inter layer Coulomb interaction |14[[T5] . 

In this letter we develop a unified approach to transport 
properties of chiral-symmetric multilayer setups shown in 
fig. [H Our main result is the universality of the full- 
counting statistics for the charge transport at the spec- 
tral degeneracy point (e = 0). We demonstrate that this 
universality is the consequence of the chiral symmetry of 




Fig. 1: Schematic view of a rectangular trilayer graphene sam- 
ple contacted by metallic leads. 

the tight-binding model of graphene multilayer. There- 
fore, it extends to arbitrary perturbations of the model, 
which preserve the chiral symmetry. Such perturbations 
include: weak crystal strains, arbitrary magnetic field, and 
interlayer hopping disorder. 

The universal charge transport is provided by the 
evanescent modes that originate at the metal-graphene 
boundaries. We find that the universal two-terminal con- 
ductivity of the setup is given by cr = AMe^ /irh, where 
M is the number of layers, e is the electron charge, and 
h is the Planck constant. The Fano factor at zero energy 
takes on the universal value F = 1/3. These results are 
insensitive to magnetic field, weak strain, and hopping dis- 
order, which can be gauged away from the full-counting 



p-1 



W.-R. Hannes and M. Titov 




Fig. 2: The sum of transmission probabilities, X^„^i7m(f?), at 
e = as a function of the transversal momentum q in the case 
of few layer graphene (M = 5) with Bernal (solid line) and 
rhombohedral (dashed line) stacking and vanishing magnetic 
field. 




Fig. 3: Quantum-mechanical (dashed lines) and averaged con- 
ductance of AB-stacked few-layer graphene as a function of 
the Fermi-energy. The layer number is one through five from 
bottom to top. The inset is a zoom into the low-energy region. 

statistics by means of a non- Abelian gauge transformation. 
Individual transmission eigenvalues are, however, less uni- 
versal and strongly depend on the stacking configuration 
(see fig.©. 

In the second part of the letter we adopt the scatter- 
ing approach to calculate the full-counting statistics of 
the charge transport away from the degeneracy point for 
£ ^ hv/L, where L is the distance between the metallic 
leads and v is the Fermi velocity. This transport regime is 
dominated by the propagating modes that give rise to the 
sample-dependent Fabri-Perot oscillations in the conduc- 
tance and noise of the ballistic setup. The gate-voltage 
dependence of these quantities is further analyzed by av- 
eraging over the Fabri-Perot oscillations. Our results for 
regular stacks are presented in figs.[5][S] For any multilayer 
above the highest band threshold we find the averaged con- 
ductance G = Me'^W\e\/h and seemingly universal value 
of the Fano factor, F = 1/8. 
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Fig. 4: Averaged Fano factor in transport through AB-stacked 
few-layer graphene. The horizontal line corresponds to F — 
1/8. The inset shows a comparison with the exact result 
(dashed curve) found from eqs. (|9l7l8p . 
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Fig. 5: Exact (dashed line) and averaged (solid line) Fano fac- 
tor for ABC-stacked trilayer graphene. The inset shows the 
transmission spectrum. 

Model. — The chiral symmetry is preserved in 
the nearest-neighbour tight-binding model of few-layer 
graphene for the regular stacking types, Bernal (AB- • • ) 
and rhombohedral (ABC---), and for all irregular ones 
excluding sequences of the type AA. 

For these models the Fermi surface becomes point-like 
at the spectrum degeneracy point (e = 0). At this energy 
no propagating states exist, hence the conductivity of the 
idealised material is expected to vanish. This is, however, 
not true in the setup of a finite size shown in fig. [T] due 
to the presence of zero-energy evanescent modes, which 
contribute to the transport in a peculiar way P^H2^ . The 
theory predicts the conductance of a ballistic monolayer 
setup G = aW/L, where the two-terminal conductivity, 
(T, approaches the universal value Ae^/nh at zero energy, 
while the Fano takes on the value F = 1/3 as in a diffusive 
system. Such pseudo-diffusive transport regime has been 
indeed observed in experiments [531 [H] with submicron 
monolayer flakes. 
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Here we extend and generalize these results for few- 
layer chiral-symmetric graphene. The model Hamiltonian 
is given by the sum of four terms H = Hq + Ha + Hs + H± , 
where Hq describes the system of AI uncoupled monolay- 
ers, Ha and Hs account for magnetic field and weak in- 
ternal strain, correspondingly, and H± takes into account 
the coupling between the layers. 

In the effective mass approximation the Hamiltonian of 
isolated monolayers. 
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(v) 



(1) 



consists of AI copies of the familiar two-dimensional Dirac 
equation with the velocity v ~ lO^m/s and the vec- 
tor of Pauli matrices cr = ((Jx,(Ty)- The notations 1 



and Ig^^ stand for the unit matrices in spin and valley 
space, respectively. The external magnetic field is de- 



{ev/c)cr ■ A ® 1 



1 



{v) 



scribed by the term Ha 
where A{r) = diag(Ai, • ■ • , Am) with A^ is the vec- 
tor potential in the m-th layer. The weak crystal strain 
(e.g. induced by ripples) can be taken into account by 
a similar term Hs = cr ■ S ® \''2^ ® Tz [15] with a strain 
field, S{r) = diag(S'i, ■ • ■ , S'm), where stands for 
the Pauli matrix in the valley space. To simplify the no- 
tations we omit outer products and let hv = ev/c = 1 in 
those expressions where it cannot cause confusion. 

The interlayer coupling is described by taking into ac- 
count the nearest-neighbour hopping both for intralayer 
and interlayer processes }26l[?7] . 
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with t±^ = 0.3-0.4 eV and s„ = cr^ = (Ta; ± lay, where the 
choice of sign depends on the stacking order. The ener- 
getically most favourable configuration is the Bernal stack- 
ing (AB) characterised by S2m = tT+, S2m+i = cr_. The 
less favourable rhombohedral configuration (ABC) corre- 
sponds to Sm = 0-+ for TO = 1, 2, ... M — 1. The model is 
readily generalised to account for spatial variations of the 
hopping parameter. 

The central property of our model is the chiral symme- 
try 

UzHa, = -H, (3) 

which holds for any spatial dependence of the vector po- 
tential A, strain field S, and the hopping parameter t±. 

It follows from the chiral symmetry ^ that the spec- 
trum determined by the "Dirac" equation, H'^ = e'^, is 
symmetric with respect to e = 0. Analytical results are 
available for regular multilayers with Ha = Hs = 0. In 



particular, the spectrum of the ABC-stacked multilayer is 
gapless and contains two branches touching at e = , 



(4) 



where £± = hv/t± = 1.6-2.2nm and fc is the wave-vector. 
The spectrum of AB-stacked multilayer is given by [5S] 



— t I C„, ± t I 



COS ■ 



M+V 



(5) 



where m = 1, 2, . . . , M. The result ^ is due to the exact 
mapping [55] onto a system of A'I/2 bilayers (for even AI) 
or onto a system of a single monolayer and {AI — l)/2 
bilayers (for odd AI). 

For both models (|4|) and ([5]) the Fermi surface at e = 
becomes point-like. We shall see that this property ex- 
tends to a general model with arbitrary spatially depen- 
dent terms Ha, Hs and H^^. 

We define the matrix Cl{r), which yields the equation 
— crVfi = Ha + Hs + H± and introduce the local gauge 
transformation [^ 



(6) 



Using the chiral symmetry ^ one can show that the 
transformation ([5]) relates the zero-energy eigenstate, \l/o, 
of the full Hamiltonian, if^'o = 0, to the eigenstate, 3>o, 
of the bare system, Hq^q = 0. Since the spectrum of Hq 
is gapless at e = and the Fermi surface is degenerate at 
e = 0, the same properties apply to the full model, H. 

In the following sections we explore the consequences 
for the transport properties of graphene multilayer. 

Transfer matrix approach. We apply the scat- 
tering approach to the transport in graphene multilayers 
for the two-terminal setup depicted in fig. [T] disregarding 
interaction effects. In particular, we calculate the full- 
counting statistics for the charge transport at zero fre- 
quency and temperature, which is determined by the cu- 
mulant generating function j30j 



J"(x) = Indet [1 - iP 



(7) 



where t is the matrix of transmission amplitudes for the 
scattering states in the leads at a given energy e. The 
conductance and the Fano factor are given by, respectively. 



G = Gq lim ^ , 



o2 77 

F^iGo/G)lun—, (8) 



where Go = 4e2//i is the conductance quantum (the factor 
of 4 takes into account the spin and the valley degeneracy). 

Following earlier works [TB] [181 US] we apply sharp 
boundary conditions at the metal-graphene interfaces at 
a; = and x = L. We also restrict ourselves to the limit 
W ^ L, where W is the width and L is the length of the 
rectangular graphene sample. Without loss of generality 
we apply periodic boundary conditions in y direction, (a 
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particular form of the boundary conditions is irrelevant 
for W ^ L), hence the momentum quantization in the 
leads, Qn ~ 2T:n/W, where q stands for the projection of 
the momentum on y axis. 

The transport properties are, then, characterised by 
the evolution matrix, T, which relates the wave-function 
Fourier components in the left and the right lead ^q(i) = 
J2q' 7qq''^q'{0)- In the casB of rigid boundary conditions 
the relation between the evolution matrix and the transfer 
matrix, A4, takes the simple form [19], Ai = C^TC with 
£■ = {ffz + ax)/V2- Using the relation 



MM^ 



W 








we rewrite the cumulant generating function ([7]) as 
1 



J" 



■ In det 



I -TV 



4ex rr^ 



(9) 



(10) 



In the next section we use this form of the generating func- 
tion to analyse the transport properties of chiral graphcnc 
multilayers at £ = 0. 

Transport at zero energy. — The transformation 
(|6]) enables us to calculate the charge transport at e = 
for a general model. Let us start, however, from a simple 
case of ballistic multilayer in the absence of magnetic and 
strain fields. Ha = Hs = 0. In this case, the gauge trans- 
formation dH]) reduces to vl* = e'^^$ with E = -axH±. 
Using this transformation we obtain the evolution matrix 



T 



7-7-T 



e-^qLp-i 



(11) 



iSL„-iS + L 



We used that the ma- 



with diag(P,P-i) = e'^^e 
trices az and S commute. The eigenvalues of the matrix 
P are parameterised by exp(— 2KmL). Since the eigen- 
values coincide with those of P~^, they appear in pairs 
with Km,m' = i|«^m|- (Thc Unpaired eigenvalue for odd 
M corresponds to Km = 0.) From equations ((TT|) . ^ one 
finds individual transmission probabilities (which are the 
eigenvalues of tt'') as 



Tniiq) = [cosh^((7- Km)L] 



1,2, ...Af. (12) 



The values of Km play the role of momentum shifts, which 
are irrelevant in the limit W ^ L, since the quantiza- 
tion of q is dense. In this limit the zero-energy generating 
function is obtained from equation (jlO[) as 



J^{x) = -M{W/ttL) [arccos(expx/2)] 



(13) 



This form of the generating function coincides with that 
of a diffusive system [21] despite the ballistic nature of 
charge transport in our model. The direct consequence of 
equation (|13[) is the universal form of the conductivity, G, 
and the Fano factor, F, at zero energy 



G = Ae'^MW/irhL, F = 1/3. 



The full-counting statistics and the results hold 
irrespective of the stacking order between the layers and 
are independent of the stacking specific momentum shifts 
Km- These results rely on the validity of the sharp bound- 
ary conditions with max Km ^ 27r/Ai?, where Xp stands 
for the Fermi wave length in thc lead. 

The momentum dependence of the transmission proba- 
bilities is less universal which is illustrated in fig. jS] For a 
multilayer with AB stacking we find 



[cmL/li_ + Vl + cUL/e±r) , (15) 



where the coefficients Cm are defined in eq. ([5]) and m = 
1,2, . . . , M . This result is consistent with thc mapping of 
the multilayer onto independent bilayers [26]. For M = 2, 
the result of eq. (jl5|) has been obtained by Snyman and 
Beenakker [T5] . 

For the ABC stacking configuration wc find the asymp- 
totic expressions 



= L 



\2m-M-l 



(16) 



(14) 



in the limit £± <^ L. The transmission resonances shown 
in fig. 121 are much better separated for the ABC multilayer 
than for thc multilayer with the Bernal stacking. The mo- 
mentum shifts Km in both cases depend logarithmically on 
the ratio L/£j_, so that the condition maxK„i <^ 2tt/Xf, 
is hard to violate. The validity of eq. ((T2|) is restricted to 
|e| ^ Eq, where Sq = t±/ {2ci){n£±/ L)'^ for AB and Eq « 
t±{TT£±/ L)^^ for ABC stacking. For ballistic graphcnc rib- 
bons with W < L, thc full counting statistics is sensitive 
to the shifts Km due to thc transversal momentum quan- 
tization. 

Remarkably, the results (|13I14|) remain valid even in 
the presence of arbitrary magnetic and strain fields. To 
justify this statement it is convenient to consider the evo- 
lution operator T in the real space representation, such 
that '^{L, y) = T^'(0, y). Using the transformation ^ we 
find 

where the magnetic and strain fields and inter-layer cou- 
pling are entering solely by means of the matrix phase 
taken at the graphene-metal boundaries. The gauge of 
thc vector potential can be chosen such that the phase 
at the boundary is y-independent. Hence the matrix ex- 
ponents in eq. (|T7|) commute and the evolution opera- 
tor in the channel space takes the form of eq. (fTTj) with 

E = {(l{L) - n{o))/L. 

Gate- voltage dependence. In the vicinity of the 
Dirac point, for e <C hv/L, the transport is entirely due 
to the evanescent modes, which are responsible for the 
pseudo-diffusive form of the full counting statistics (fT3)) 
at e = 0. Away from the Dirac point, the transport is 
less universal and depends on a number of details. In 
this section we restrict ourselves to ballistic models with 
Ha=Hs = 0. 
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We calculate the transmission probabilities Tm{q) and 
perform the numerical integration over the transverse mo- 
mentum q (assuming 2tt/Xf 3> |£|) to find the conduc- 
tance and the Fano factor in the limit W ^ L. In the 
case of AB stacking the probabilities Tm{q) can be found 
analytically and the resulting energy dependence of the 
conductance is shown in fig. [3] with dashed lines. The 
Fano factor for the trilayer with AB and ABC stacking is 
depicted with the dashed lines in fig. |4] (inset) and fig. [5l 
correspondingly. At £ = the figures confirm the results 
of eqs. (fnm)) . 

For energies exceeding the ballistic Thouless energy, 
hv/L, the transport is dominated by propagating modes, 
which give rise to the sample-specific Fabry-Perot oscil- 
lations in conductance and noise. In order to get exper- 
imentally relevant results we perform the averaging over 
these oscillations treating the propagating phases as ran- 
dom quantities, which are uniformly distributed in the in- 
terval (0, 27r) This type of averaging corresponds 
to a quasiclassical approximation that respects the con- 
servation of the transversal momentum, q, in the sample. 

We introduce the individual scattering matrices for the 
left (L) and right (R) sample-lead interfaces 



L/R 



PL/R Tl/r, 



(18) 



which relate the scattering state amplitudes in the leads 
with those in the sample. We assume that the number 
of open channels in the sample equals Mq < M for a 
given energy, e, and transversal momentum, q. In the 
planar geometry of fig. [T] the S'-matriccs in eq. (|18p are 
readily calculated by matching the scattering states in the 
corresponding lead and in the sample. 

As the result the total transmission matrix from the left 
to the right lead can be written as 



i = TR (1 - -To Pl To Pr) ^ To TL, 



(19) 



where fo = diag(e"^^ , • • • , e'*"o ) is parameterised by the 
propagating phases (j)m = k„iL accumulated in the free 
propagation inside the sample. The number of propa- 
gating channels, Mq, is determined from the requirement 
Im(A:,„) = 0, where the longitudinal momentum km is 
found from the dispersion relation s — e{k,q) for a given 
energy s and transversal momentum q. 

An equivalent way to calculate the averaged conduc- 
tance is formulated in terms of the classical transmission 



and reflection probabilities, TRYL^ 



^, and 



Ri 



R(L),' 



Ipr(l) 



for the transport through the 



sample-graphene interfaces. The classical probabilities to 
pass through an entire sample are, then, organised in the 
following matrix T = Tr(1 — RlRr)~^Tl, for each value 
of q. The quantity, T, ignores the phase-coherence in 
the assumption that the transversal momentum is con- 
served inside the sample. The averaged conductance, 
G = Go Eg E„™ T„„ coincides with G = Go E,(Trtt^), 



where the brackets stay for the averaging over the propa- 
gating phases 0™ in Eq. p^ . 

A special case of a single propagating channel, Mq ~ 1, 
per transversal momentum, g, is naturally realised in a 
monolayer and multilayers with a single ungapped band 
for energies below the lowest band threshold. The example 
of the latter is the ABC-stacked multilayer. For equivalent 
sample-lead junctions, the cumulant generating function 
([7]) can be expressed using eq. (jl9|) as 



T^ex + 4(1 - T) sin^ 
T2 + 4(1-7') sin% 



(20) 



where T = f f ^ is the g-dependent transmission probabil- 
ity of the sample-lead interface. The averaging over the 
propagating phase, 0, can be performed analytically with 
the result 



J"(x) = 2 ^ In Tc^^^ + y^4(l -T) + T^cx 



(21) 



With the help of eq. ([U we obtain for the averaged con- 
ductance and the Fano factor 



G = GoJ2 



T 



2-T 



F 



Go ^ 2T(1 - f) 

G ^ (2 - r)3 ■ 



(22) 



Let us digest eqs. (j21l22[) in the case of monolayer 
graphene. If the leads in the setup of fig. Q are modelled 
by highly doped graphene one finds the interface transmis- 
sion probability as T{q) = 1 — tan^ 9/2, where q = \e\ sin 6' 
and 6 e (— 7r/2, tt/2) is the angle of incidence. For W ^ L 
wc replace the summation in eqs. (j22p with the integration 
over q in the interval (— e, e) and reproduce the asymptotic 
results EOj G = GoW\s\/4: and F = 1/8. 

Similarly, for the bilayer we find from eqs. (|21l22p the 
averaged conductance below the band threshold. 



G = GoW ( Y 



t 



2 7-'? 



e <ti 



(23) 



The result 



where 7=1 + 2|£|/i^ and 77 ~ yl+T — 7 
for the Fano factor in the bilayer setup takes the form 



GqW f\e\ ^ 4c^+C2/hhj) 
G [32^ ^ 512|e|3 



|e|<U, (24) 



where ci = 3 — 77(1+7) — 7'' and C2 = 3 + 7(7 + 2) (167'' + 
1273 - 1772 + IO7- 6). 

If several propagating channels per the value of q open 
up in the sample, the averaging procedure is compli- 
cated and has to be carried out numerically. Still, for 
AB-stacked multilayers the analysis is simplified by us- 
ing Koshino-Ando mapping [26j to an effective bilayer- 
monolayer system. For energies below the first band 
threshold, the conductance and noise can be constructed 
from the available analytical results for mono- and bi-layer 
graphene using the effective bilayer coupling constants. 
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tf,m = '^c,nt±, where m = 1, 2, . . . , Int[A//2]. For higher 
energies, one finds at most two propagating channels per q 
so that the numerical implementation of the averaging pro- 
cedure is straightforward. The band thresholds are seen 
as kinks in the energy dependence (transmission spectra) 
of conductance and noise in figs. |3] and |31 At e = wc for- 
mally find G = and (for AB-stackcd multilayer) F = 1/2 
for even M and F = (6M - 5)/(12M - 4) for odd mmr- 
bcr of layers, M . These results ignore the contribution of 
evanescent modes. 

The interlayer coupling becomes irrelevant for the trans- 
port at energies far above all band thresholds. Universal 
asymptotic results, G = MGoW\e\/4: and F = 1/8, are 
obtained in this limit for any combination of AB and ABC 
stacking. This universality is due to the linear dispersion, 
e Cmt± ± fi,i;|fe|, of all spectral branches of our effective 
few-layer model at high energies. 

One can modify the decomposition (jl9p to account for 
evanescent modes, which correspond to imaginary values 
of the propagating phases, 4>m = kmL, and use an appro- 
priate analytic continuation of the matrices ^l/r, which 
become non-unitary. We employ this approach to plot 
the exact conductance and Fano factor for ABC-stacked 
graphenc in fig. [SJ 

Conclusions. — A class of Hamiltonians describing 
few-layer graphenc with spatially dependent strain and 
magnetic field obeys a chiral symmetry. This symmetry is 
responsible for the universal pseudo-diffusive charge trans- 
port via evanescent modes. The universality is understood 
on the basis of a non-Abelian gauge transformation at zero 
energy. For finite energies the full counting statistics of 
few-layer ballistic samples is calculated using the scatter- 
ing approach. 
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